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Abstract 

In this paper we address some modelling issues related to biological growth. Our treatment 
is based on a recently-proposed, general formulation for growth within the context of Mixture 
Theory (Journal of the Mechanics and Physics of Solids, 52, 2004, 1595-1625). We aim to en- 
hance this treatment by making it more appropriate for the biophysics of growth in porous soft 
tissue, specifically tendon. This involves several modifications to the mathematical formulation 
to represent the reactions, transport and mechanics, and their interactions. We also reformulate 
the governing differential equations for reaction-transport to represent the incompressibility con- 
straint on the fluid phase of the tissue. This revision enables a straightforward implementation 
of numerical stabilisation for the hyperbolic, or advection-dominated, limit. A finite element 
implementation employing an operator splitting scheme is used to solve the coupled, non-linear 
partial differential equations that arise from the theory. Motivated by our experimental model, 
an in vitro scaffold- free engineered tendon formed by self-assembly of tendon fibroblasts ( Tissue 
Engineering, 10, 2004, 755-761), we solve several numerical examples demonstrating biophysical 
aspects of tissue growth, and the improved numerical performance of the models. 



1 Introduction 



Growth involves the addition or depletion of mass in biological tissue. Growth occurs in combi- 
nation with remodelling, which is a change in microstructure, and possibly with morphogenesis, 
which is a change in form in the embryonic state. The physics of these processes are quite distinct , 
and for modelling purposes can, and must, be separated. Our prev ious work (Garikipati et all 
2004 ), upon which we now seek t o build, drew in some measu re from Cowin and Hegedus ( 19761 ): 
Epstein and Maugin ( 2000l ). and Taber and Humphrey ( 2001 ). and was focused upon a compre- 
hensive account of the coupling between transport and mechanics. The origins of this coupling 
were traced to the balance equations, kinematics and constitutive relations. A major contribution 
of that work was the identification and discussion of several driving forces for transport that are 
thermodynamically-consistent, in the sense that specification of these relations does not violate the 
Clausius-Duhem dissipation inequality. 

There have been a number of significant papers on biological growth and remodelling in the last 
7-8 years of which we touch upon some, whose approach e s are either similar to ours in some respects 
or differ in important ways. Humphrey and Raiagopal ( 20021 ) provided a mathematical treatment 
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of adaptation in a tissue, which includes growth and remodelling in the sense of this paper. The 
authors identified adaptation as perhaps the most important mechanical characteristic of biological 
tissue. They introduced the notion of evolving natural configurations to model the state of material 
deposited at different instants in time. The treatment of the growth part of the deformation 
gradient in this paper bears some resemblance to this idea, although a detailed development has 
not been pursued here. The focus, instead, is on detailing some aspects of the problem that 
derive from treatment of the tissue as a porous medium, or as a mixture of interacting species. 
Preziosi and Farina ( 20021 ) developed an extension to the classical Darcy's Law to incorporate mass 
exchanges between reacting species. This consideration is relevant to growth problems; however, in 
our opinion, these issues were subsumed in iGarikipati et al.1 (|2004l ). upon which this paper is based. 
Many of the ideas employed here are applicable to tumour growth problems; however, due to our 
current focu s on tendon, we do no t include phenomena such as angiogenesis and cell migration (see 



for example iBrewardet aH B. The changes in concentration that occur with growth tend to 



cause swelling or contraction of the tissue. This phenomenon has been accounted for pre viously by 



us in fields unrelated to Bio l ogy, u sing the idea of thermal expansion. See, for example, I Rao et al. 



(l2000h and lOarikipatT et al.l (|200l l). which, too, are probably no t the first instance s of th is idea. In 
the literature on biological growth this connection was made by lKlisch and Hoger (|2003h . 

In the present paper, we seek to restrict the range of physically- admissible models in order to 
gain greater physiological relevance for modelling growth in soft biological tissue. We also include 
one improvement in the mathematical/numerical treatment: The advection-diffusion equations for 
mass transport require numerical stabilisation in the advection-dominated regime (the hyperbolic 
limit). We draw upon the enforcement of the incompressibility limit for the fluid phase to facilitate 
this development. Below, we briefly introduce each aspect that we have considered, but postpone 
details until relevant sections in the paper. 

• For a tissue undergoing finite strain, the transport equations can be formulated, mathemat- 
ically, in terms of concentrations with respect to either the reference or current (deformed) 
configuration. However, the physics of fluid-tissue interactions and the imposition of relevant 
boundary conditions is best understood and represented in the current configuration. 

• The state of saturation is crucial in determining whether the tissue swells or shrinks with 
infusion/expulsion of fluid. This aspect has been introduced into the formulation. 

• The fluid phase, whether slightly compressible or incompressible, can dev elop compressiv e 
stress without bound. However, it can develop at most a small tensile stress ( Brennen . 19951 ). 
having implications for the stiffness of the tissue in tension as against compression. Although 
this also has implications for void formation through cavitation, the ambient pressure in 
the tissue under normal physiological conditions ensures that this manifests itself only as a 
reduction in compressive pressure. 



Whe n modelling transport, it is common to assume Fickean diffusion (IKuhl and Steinmann . 
20031 ). This implies the existence of a mixing entropy due to the configurations available to 
molecules of the diffusing species at fixed values of the macroscopic concentration. The state 
of fluid saturation directly influences its mixing entropy. 

If fluid saturation is maintained, void formation in the pores is disallowed even under an in- 
crease in the pores' volume. This has implications for the fluid exchanges between a deforming 
tissue and a fluid bath with which it is in contact. 
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Recognising the incompressibility of the fluid phas e, it is comm on to treat soft biological tissue 
as either incompressible or nearly- incompressible ( Fung . 19931 ), At the scale of the pores (the 
microscopic scale in this case), however, a distinction exists in that the fluid is exactly (or 
nearly) incompressible, while the porous solid network is not obviously incompressible. 



• In Garikipati et al. ( 20041 ). the acceleration of the solid phase was included as a driving force 
in the constitutive relation for the flux of other phases. However, acceleration is not frame- 
invariant and its use in constitutive relations is inappropriate. 

• Chemical solutes in the extra-cellular fluid are advected by the fluid velocity and additionally 
undergo transport under a chemical potential gradient relative to the fluid. In the hyperbolic 
limit, where advectio n dominates, spatial instabilities emerge in nu merical solutions of these 
transport equations ( Brooks and Hughes . 19821 ; Hughes et al. . 19871 ) . Numerical stabilisation 
of the equations is intimately tied to the mathematical representation of fluid incompressibility. 

• The modelling of solid-fluid mechanical coupling carries strong implications for the stiffness 
of tissue response, the nature of fluid transport, and since nutrients are dissolved in the fluid, 
ultimately for growth. We present upper and lower bounds for this problem and computations 
of coupled boundary value problems with these bounds. 

These issues are treated in detail in relevant sections of the paper, which is laid out as follows: 
Balance equations and kinematics are discussed in Section [2j constitutive relations for reactions, 
transport and mechanics in Section (3J and numerical examples are presented in Section HI Conclu- 
sions are drawn in Section [5j 



2 Balance equations and kinematics of growth 

In this section, the coupled, continuum balance equations governing the behaviour of growing tis- 
sue are summarised and specialised as outlined in Section [TJ For detai led continuum mechan ical 
arguments underlying the equations, the interested reader is directed to Garikipati et al. ( 20041 ) . 

The tissue of interest is an open subset of M 3 with a piecewise smooth boundary. At a reference 
placement of the tissue, f^o, points in the tissue are identified by their reference positions, X 6 f^o- 
The motion of the tissue is a sufficiently smooth bijective map ip : Qq x [0, T] —> M 3 , where 
Oo := U <9f2o! d£lo being the boundary of Qo- At a typical time t £ [0, T], <p(X,t) maps 
a point X to its current position, x. In its current configuration, the tissue occupies a region 
Vtt = fti^o)- These details are depicted in Figured) The deformation gradient F := d(p/dX is 
the tangent map of (p. 

N ■ M' 




Figure 1: The tissue as a continuous medium with growing and diffusing species. 
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The tissue consists of numerous species, of which the following groupings are of importance for 
the models: A solid species, consisting of solid collagen fibrils and ce//sQ denoted by c, an extra- 
cellular fluid species denoted by f and consisting primarily of water, and solute species, consisting 
of precursors to reactions, byproducts, nutrients, and other regulatory chemicals. A generic solute 
will be denoted by s. In what follows, an arbitrary species will be denoted by t, where i = c, f , s. 

The fundamental quantities of interest are mass concentrations, p^X^t). These are the masses 

of each species per unit system volume in Qq. Formally, these quantities can also be thought of 

in terms of the maps p L : x [0, T] — > R, upon which the formulation imposes some smoothness 

requirements. By definition, the total material density of the tissue at a point is a sum of these 

concentrations over all species ^ p L = pq. Other than the solid species, c, all phases have mass 
_ l 

fluxes, M'B These are mass flow rates per unit cross-sectional area in the reference configuration 
defined relative to the solid phase. The species have mass sources (or sinks), IT. 



2.1 Balance of mass for an open system 



As a result of mass transport (via the flux terms) and inter-conversion of species (via the source/sink 
terms) introduced above, the concentrations, p^, change with time. In local form, the balance of 
mass for an arbitrary species in the reference configuration is 



IT -DIV[M l ], W, 



(1) 



M 

dt 

recalling that, in particular, M c = 0. Here, DIV[«] is the divergence operator in the reference 
configuration. The functional forms of IT 1 are abstractions of the underlying biochemistry, physio- 
logically relevant examples of which are discussed in Section \3. 41 and the fluxes, M L , are determined 
from the thermodynamically-motivated constitutive relations described in Section 13.31 

The behaviour of the entire system can be determined by summing Equation ([1]) over all species 
i. Additionally, sources and sinks satisfy the relation 



£IT 



0, 



(2) 



which is consistent ( Garikipati et al. . 20041) wi th the Law of Mass Action for reaction rates and 
with Mixture Theory (jTruesdell and Nolll . Il965l ) . 



2.1.1 The role of mass balance in the current configuration 

In order to proceed, we must first introduce the central kinematic assumption underlying the for- 
mulation: We assume that the pore structure deforms with the collagenous phase. Therefore, the 
deformation gradient, F, is common to c and the fluid-filled pore spaces. Furthermore, in what 
follows, we will treat the fluid as ideal and nearly-incompressible, i.e. as elastic (Section 13. 2p . This 
combination of kinematic and constitutive assumptions to be elaborated upon, implies that the 
stress in the fluid phase is determined by the elastic part of F (see Sections 12.21 and I3.2() . For clarity 
we denote it as F c . Importantly, the pore-filling fluid under stress can also undergo transport 
relative to the pore network; i.e., relative to the collagenous phase. This is the fluid flux, denoted 

At this point, we do not distinguish the solid species further. This is a good approximation to th e physiological 
settin g for tendons, which are relatively acellular and whose dry mass consists of up to 75% collagen (|Nordin et all . 
l200lft . 

2 Currently, we do not consider certain physiological processes, such as the migration of fibroblasts within the 
extra-cellular matrix during wound healing, which may otherwise be modelled as mass transport. 



4 



rt i n 



Figure 2: If the pore structure at the boundary deforms with the tissue and this boundary is in 
contact with a fluid bath, the fluid concentration with respect to the current configuration, i.e., //, 
remains constant. 



by M in the reference configuration. Note that the specification of constitutive relations for the 
flux is still open at this point in the discussion. At the outset, we preclude stress in any of the 
solute species, s. Only the solid collagen and fluid bear stress. 

Although the initial/boundary value problem of mass transport can be consistently posed in the 
reference configuration, the evolving current configuration, 0^, is of greater interest from a physical 
standpoint for growth problems. It follows from the discussion in the preceding paragraph that the 
shape and size of pores in Of is determined by F. Therefore, at the boundary, the fluid concentration 
with respect to Qt remains constant if the boundary is in contact with a fluid bath. Accordingly, this 
is the appropriate Dirichlet boundary condition to impose under normal physiological conditions. 
This is shown in an idealised manner in Figure [2J 

In the interest of applying boundary conditions (either specification of species flux or concen- 
tration) that are physically meaningful, we use the local form of the balance of mass in the current 
configuration, 

^ = tt l - div[m l ] - p L dw[v] , V l, (3) 

where p L (x , t) , tt l (x , t) , and m i (x,t) are the current mass concentration, source and mass flux of 
species i respectively and v(x, t) is the velocity of the solid phase. They are related to corresponding 
reference quantities as p L = (det (F))" 1 po> = (det (F)) -1 IP and m l = (det (F)y 1 FM L . The 
spatial divergence operator is div [•], and the left hand-side in Equation (J3j) is the material time 
derivative relative to the solid, which may be written explicitly as -§i\x, implying that the reference 
position of the solid collagenous skeleton is held fixed. 



2.2 The kinematics of growth induced by changes in concentration 



Local volumetric changes are associated with changes in the concentrations of the solid collagen 
and fluid, i = c, f . If the material of the solid collagen or fluid remains stress free, it swells 
with an increase in concentration (mass of the species per unit system volume), and shrinks as its 
concentration decreases. This leads to the notion of a growth deformation gradient. One aspect of the 
coupling between mass transport and mechanics stems from this phenomenon. In the setting of finite 
strain kinematics, the total deformation gradient, F, is decomposed into the growth component of 
the solid collagen, F e , a geometrically-necessitated elastic component accompanying growth, F 

— e c c — e c ~ e c 

and an additional elastic component due to external stress, F . Later, we will wr ite F e = F F . 



This split is analogous to the classical decomposition of multiplicative plasticity (jLed . Il969f) and is 



similar to the approach followed in existing literature on biolog ical growth (see for e.g. iKlisch et al 



200 ll ; iTaber and Humphrevl l200ll ; lAmbrosi and Mollical . [2002]) . As explained in Section 12.1.11 we 



assume that the fluid-filled pores also deform with F, and that a component, F c , of this total 
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Figure 3: The kinematics of growth, which holds for i = c, f . 

deformation gradient tensor, determines the fluid stress. We also assume a fluid growth component, 

F g , which we elaborate below, and that F° F g = F. As with the solid collagen we admit 

f e f ~ e f 

F c = F F , the sub-components carrying the same interpretation as for the solid collagen. 
However, we do not explicitly use this last decomposition. 

The elastic-growth decomposition is visualised in Figure Assuming that the volume changes 
associated with growth described above are isotropic, a simple form for the growth part of the 
deformation gradient tensor is 

(,:;:) 1 ^ (4) 

where pfa (X) is the reference concentration at the initial time, and 1 is the second-order isotropic 
tensorH In the state, F = F g , the species would be stress free. The kinematics being local, the 
action of F g alone can result in incompatibility, which is eliminated by the geometrically-necessary 

elastic deformation F , which causes an internal, self-equilibrated stress. The component F is 
associated with the external stress. 



2.2.1 Saturation and tissue swelling 

The degree of saturation of the solid phase plays a fundamental role in determining whether the 
tissue responds to an infusion (expulsion) of fluid by swelling (shrinking). In particular, the isotropic 
swelling law defined by Equation @ has to be generalised to treat the case in which the solid phase 
is not saturated by fluid. 

Figure 0] schematically depicts two possible scenarios. If the tissue is unsaturated in its current 
configuration, as in A, then, on a microscopic scale, it contains unfilled voids. It is thus capable 
of allowing an influx of fluid, which tends to increase its degree of saturation until fully saturated, 
as in B. This increase does not cause swelling of the tissue in the local stress-free state, as there 
is free volume for incoming fluid to occupy. However, once the tissue is saturated in the current 

3 This choice is only the simplest possible. Given the highly directional micro-structure and mechanical properties 
of many tissues, it seems likely that anisotropic growth is actually more common. Wolff's Law for bone growth is 
one example. This is a topic of ongoing investigation, and one that we will report on in greater detail in a future 
communication. 
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Figure 4: Unsaturated tissue in the current configuration (A) allows influx of fluid without swelling 
until it is completely saturated (B). Initially saturated tissue (B), in general, swells with influx of 
fluid (C). 

configuration, an increase in the fluid content causes swelling in the stress-free state, as depicted 
in C, since there is no free volume for the entering fluid to occupy. It is this second case that is 
modelled by (|4|). It is worth emphasizing that this argument holds for F s , which is the local stress- 
free state of deformation of the fluid-containing pores at a point. The actual deformation gradient, 

f f f 
F = F e F g , also depends on the the elastic part, F e , which is determined by the constitutive 

response of the fluid. Under stress, an incompressible fluid will have det.F c = 1 and therefore a 

f 

fluid-saturated tissue will swell with fluid influx, det.F = det.F g > 1. A compressible fluid may 
have detF e < 1 allowing det.F < 1 even with det.F g > 1. Even in this case, however, in the 
stress-free state there will be swelling. 

Therefore, for the fluid phase, the isotropic swelling law can be extended to the unsaturated 
case by introducing a degree of saturation, v\ defined in the current configuration, f^. We have 
v L = p L /p L , where p l is the intrinsic density in Qt and is given by p L = p^/detF. Note that the 
intrinsic reference density, /5q, is a material property. Upon solution of the mass balance equation 
([3|) for p\ the species volume fractions, v\ can therefore be computed in a straightforward fashion. 
The sum of these volume fractions is our required measure of saturation defined in f^. Also, 
recognizing that for the dilute solutions obtained with physiologically-relevant solute concentrations, 
the saturation condition is very well approximated by v f + v s = 1, we proceed to redefine the fluid 
growth-induced component of the pore deformation gradient tensor as follows: 

1, otherwise. 

In ([5]) pg sat is the reference concentration value at which the tissue attains saturation in the current 
configuration. 

f ~ 

With this redefinition of F g it is implicit that v + v s > 1 is non-physical. Saturation holds in 
the sense that u f + v s = 1, and it actually allows s ^2 li v L > 1 if the sum is over all species. It has 
been common in the soft tissue literature to assume that, under normal physiological conditions, 
soft tissues are fully saturated by the fluid and Equation is appropriate for l = f. However, 
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this treatment of saturation and swelling induced by the fluid phase is necessary background for 
Section 13.2.11 where we discuss the response of the fluid phase under tension. This treatment also 
holds relevance for partial drying, which ex vivo or in vitro tissue may be subject to under certain 
laboratory conditions, and is central to the mechanics of drained porous media other than biological 
tissue, most prominently, soils. 



2.3 Balance of momenta 

In soft tissues, the species production rate and flux that appear on the right hand-side in Equa- 
tions (HJ) and ([3j) , are strongly dependent on the local state of stress. To correctly model this 
coupling, the balance of linear momentum should be solved to determine the local state of strain 
and stress. 

The deformation of the tissue is characterised by the map <f(X, t). Recognising that, in tendons, 
the solid collagen fibrils and fibroblasts do not undergo mass transport, the material velocity of this 
species, V = d<p/dt, is used as the primitive variable for mechanics. Each remaining species can 
undergo mass transport relative to the solid collagen. For this purpose, it is useful to define the 
material velocity of a species i relative to the solid skeleton as: V b = (1/ Pq)FM l . Thus, the total 
material velocity of a species i is V + V 1 . 

The total first Piola-Kirchhoff stress tensor, P, is the sum of the partial stresses P L (borne by 
a species i) over all the species present. Recognizing that solutes in low concentrations, and do not 
bear appreciable stress, the partial stresses and momentum balance equation are defined only for 
the solid collagen and fluid phases. With the introduction of these quantities, the balance of linear 
momentum in local form over Qq for solid collagen and fluid is, 

Po§- t (V + n =P' (9 + Q b ) + DIV[P l ] 

- (GRAD [V + V b ]) M L , t = c,f 

where g is the body force per unit mass, and q L is an interaction term denoting the force per unit 
mass exerted upon i by all other species present. The final term with the (reference) gradient denotes 
the contribution of the flux to the balance of momentum. In practise, the relative magnitude of the 
fluid mobility (and hence flux) is small, so the final term on the right hand side of Equation (JHJ) 
is negligible, resulting in a more classical form of the balance of momentum. Furthermore, in 
the absence of significant acceleration of the tissue during growth, the left hand-side can also be 
neglected, reducing (|6|) to the quasi-static balance of linear momentum. 

The balance of momentum of the entire tissue is obtained by summing Equation §6§ over i = c, f . 
Additionally, recognising that the rate of change of momentum of the entire tissue is affected only 
by external agents and is independent of internal interactions, the following relation arises. 

f 

Y J (PW + ^V L ) = 0. (7) 



This i s also consistent with Classical Mixture Theory (jTruesdell and Nolll . ll965l ). See lGarikipati et al 



(|2004l ) for further details on balance of linear momentum, and the formulation of balance of angular 
momentum. We only note here that the latter principle leads to a sym metric partial Cauchy stress , 
cr L for each species in contrast with the unsymmetric Cauchy stress of Epstein and Maugin ( 2000l ) . 
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3 Constitutive framework and specific models 



As is customary in field theories of continuum physics, the Clausius-Duhem inequality is obtained 
by multiplying the Entropy Inequality (the Second Law of Thermodynamics) by the temperature 
field, 9, and subtracting it from the Balance of Energy (the First Law of Thermodynamics). We 
assume the Helmholtz free energy per unit mass of species t to have the formQ ifi 1, = r ^ i {F e , 6, p^). 
Substituting this in the Clausius-Duhem inequality results in a form of this inequality that the 
specified constitutive relations must not violate. O nly the valid con s tituti ve laws relevant to the 
examples that follow are listed here. For details, see iGarikipati et al.1 ((20041 ) . 



3.1 An anisotropic network model based on entropic elasticity 

The partial first Piola-Kirchhoff stress of collagen, modelled as a hyperelastic material, is P c = 

Pfjdifj j 'dF c . Recall that F e = FF g is the elastic part, and F g is the growth part, respec- 
tively, of the deformation gradient, of collagen. Following Equation (J3J), if we were considering 
unidirectional growth of collagen along a unit vector e, we would have -F gC 

denoting the initial concentration of collagen at the point. 

The mechanical response of tendons in tension is determined primarily by their dominant struc- 
tural component: highly oriented fibrils of collagen. In our preliminary formulation, the strain 
energy density for collagen has been obtained from hierarchical multi-scale considerations based 
upon an entropic elasticity-based worm- like chain (WLC) model ((Kratkv and Porodl . Il94fll ). The 
WLC model has been widely used for long chain single molecule s, most prominently for DNA 



-p-e (g) e, with pZ 



( Marko and Siggia . 1995 ; Rief et al. . 1997 ; Bustamante et al. . 20031 ). and recently for the collagen 



monomer (jSun et all 120021 ) . The central parameters of this model are the chain's contour length, L, 
and persistence length, A. The latter is a measure of its stiffness and given by A = x/kQ, where \ is 



the be nding rigidity, k is Boltzmann's constant and 6 is the temperature. See lLandau and Lifshitz 
(|195ll ) for general formulation of statistical mechanics models of long chain molecules. 

To model a collagen network struc ture, the WLC mode l has been embedded as a single con- 
stituent chain of an eight-chain model (jBischoff et all l2002aH bh . depicted in Figure [5j Homogenisa- 
tion via averaging then leads to a continuum Helmholtz free energy function, ^ C H 



pW{f°\pI) 



Nk6 f r 2 



L 



AA \2L 4(1 -r/L) 4 
1)+ 7 1: (C cC -1 



7 c" 2 ' 3 




(8) 



4 A /2L/A yl L 4(1 - y/2A/L) 4^ 



z = io g a; a c 2 a§ 



z. 



Here, ./V is the density of chains, and a, b and c are lengths of the unit cell sides aligned with 
the principal stretch directions. The material model is isotropic only if a = b = c. 



4 Conceivably, the mass-specific Helmholtz free energy of one species could be a function of the concentration of 
other species. Ion concentrations, for instance, can determine the state of osmotic tension of certain soft tissues. 
Therefore, this choice represents a constitutive restriction. 

5 Under the isothermal conditions assumed here, ?/> c is independent of 8 in the strain energy. Accordingly, we have 
the parametrisation ifj c — ^(F ,pg). 
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a 



A 



Figure 5: The eight-chain model incorporating worm-like chains. 
The elastic stretches along the unit cell axes are, respectively, denoted by Af , A| and A3, C e = 

c T c 

F e F e is the elastic right Cauchy-Green tensor of collagen. The factors 7 and (3 control the bulk 

compressibility of the model. The end to end chain length is given by r = }^\J a 2 Af 2 + 6 2 A| 2 + c 2 A| 2 , 

where Af = y Nj ■ C° c Nj, and Nj, I = 1,2,3 are the unit vectors along the three unit cell axes, 
respectively. In our numerical simulations that appear below in Section EH the numerical values 



used for the parameters introduced in ([8]) are based on those in iKuhl et allfeoOSh 



3.2 A nearly incompressible ideal fluid 

In this preliminary work, the fluid phase is treated as nearly incompressible and ideal, i.e., inviscid. 
The partial Cauchy stress in the fluid is 



cr 



f - det(F c T^F 6 = h(//)l, (9) 



where a large value of h'(p ) ensures near-incompressibility. 
3.2.1 Response of the fluid in tension; cavitation 

The response of the ideal fluid, as defined by Equation ([9]), does not distinguish between tension 
and compression, i.e., whether det(-F c ) ^ 1. Being (nearly) incompressible, the fluid can develop 
compressive hydrostatic stress without bound — a case th at is modelled accurately. However, the 
fluid can develop at most a small tensile hydrostatic stress (IBrennen , and the tensile stiffness 



6 Where, we are referring to the fluid being subject to net tension, not a reduction in fluid compressive stress from 
reference ambient pressure. 
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is mainly from the collagen phase. This is not accurately represented by ©, which models a 
symmetric response in tension and compression. 

Here, we preclude all tensile load carrying by the fluid by limiting det(F c ) < 1. We first 
introduce an additional component to the relation between deformation of the pore space, given by 
F, the fluid stress-determining tensor, F c and the growth tensor for the fluid, F g . Consider the 
cavitation (void forming) tensor, F v , defined by 

P^pg'pv = F ( 10 ) 

We restrict the formulation to include only saturated current configurations at t = 0. Following 
Section 12.2.11 we have v { + v c = 1 at t = 0, the saturation condition in fit when solutes are at 

f f 

low concentrations. At times t > Equation © holds for F g . If det [F(F g < 1 we set 
F e = F(F g ) and F v = 1 for no cavitation. Otherwise, since det[F(F g J -1 ] > 1, we specify 
F ci = det[F(F g )~ 1 ]~ 1 / 3 F(F gi )~ 1 thus restricting F e to be unimodular and allow cavitation by 
writing F v = F(F C F g ) . These conditional relations are summarized as 

{F(F gf ) _1 , F v = 1, det[F(F§V 1 ] < 1 

det[F(FS f )- 1 ]" 1/3 ^(-P 1g V 1 , ( n ) 
F v = F(F cf F g V 1 otherwise. 



3.3 Constitutive relations for fluxes 

From lGarikipati et al.1 (j2004l ). the constitutive relation for the flux of extra-cellular fluid relative to 
collagen in the reference configuration takes the following form, 



M 



D 1 



p f F T g + F T DIV 



PoGRAD/i 



(12) 



where D i is the positive semi-definite mobility of the fluid, and isothermal conditions are assumed in 
order to approxima te the physiological ones. Experimentally determ ined transport co efficients (e.g. 
for mouse tail skin ( Swartz et al. . 19991 ) and rabbit Achilles tendons ( Han et al. . 2000l )) are used for 
the fluid mobility values. The terms in the parenthesis on the right hand-side of Equation (|12p sum 
to give the total driving force for transport. The first term is the contribution due to gravitational 
acceleration. In order to maintain physiological relevance, this term has been neglected in the 
following treatment. The second term arises from stress divergence; for an ideal fluid, it reduces to 
a pressure gradient, thereby specifying that the fluid moves down a compressive pressure gradient, 



which is Darcy's Law. The third term is the gradient of the chemical potential, /r 



9rf, 



where e f is the mass-specific internal energy, 9 is temperature and rf is the mass-specific entropy. 
The entropy gradient included in this term results in classical Fickean diffusion if only mixing 
entropy exists, as discussed in the foll owing; section. For a d etailed derivation and discussion of 
Equation (fT2|) , the reader is directed to iGarikipati et al.l (|2004l ) . 



3.3.1 Saturation and Fickean diffusion of the fluid 

As depicted in Figure [6l only when pores are unsaturated are there multiple configurations available 
to the fluid molecules at a fixed fluid concentration. This leads to a non-zero mixing entropy. In 
contrast, if saturated, there is a single available configuration (degeneracy), resulting in zero mixing 
entropy. Consequently, Fickean diffusion, which arises from the gradient of mixing entropy can exist 
only in the unsaturated case. However, even a saturated pore structure can demonstrate concen- 
tration gradient-dependent mass transport phenomenologically: The fluid stress depends on fluid 
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Figure 6: Depicted at a microscopic scale, only unsaturated tissues A and B can undergo Fickean 
diffusion of the fluid. C is saturated. 

concentration (see Equation ©), and fluid stress gradient-driven flux appears as a concentration 
gradient-driven flux. 

The saturation dependence of Fickean diffusion is modelled by using the measure of saturation 
introduced in Section 12.2.11 We rewrite the chemical potential as 



M f = e f - 9r] { , 

rf -> 0, as 5 f + v c -> 1. 



(13) 



It is again important to note that under physiological conditions, soft tissues are fully saturated by 
fluid, and it is appropriate to set // = e f . 



3.3.2 Transport of solute species 

The dissolved solute species, denoted by s, undergo long range transport primarily by being advected 
by the fluid. In addition to this, they undergo diffusive transportrelative to the fluid. This motivates 
an additional velocity split of the form V s = V s + V i , where V s denotes the velocity ofjhe solute 
relative to the fluid. The constitutive relation for the corresponding flux, denoted by M s , has the 
following form, similar to Equation (|12p defined for the fluid flux. 

M s = D s (-PqGRAD [e s - Brf\) , (14) 

where D s is the positive semi-definite mobility of the solute relative to the fluid, and again, isother- 
mal conditions are assumed to approximate the physiological ones. Following Section T2.3I there are 
no stress-dependent contributions to M s . 



3.3.3 Frame invariance and the contribution from acceleration 

In our earlier treatment ( Garikipati et al. . 20041 ). the constitutive relation for the fluid flux had 



a driving force contribution arising from the acceleration of the solid phase, — PqF t ^-. This 
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term, being motivated by the reduced dissipation inequality, does not violate the Second Law and 
supports an intuitive understanding that the acceleration of the solid skeleton in one direction must 
result in an inertial driving force on the fluid in the opposite direction. However, as defined, this 
acceleration is obtained by the time differentiation of kinematic quantities^ and does not transform 
in a frame-indifferent manner. Unlike the superficially similar term arising from the gravity vectorH 
the acceleration term presents an improper dependence on the frame of the observer. Thus, its use 
in constitutive relations is inappropriate, and the term has been dropped in Equation (|12|) . 



3.3.4 Incompressible fluid in a porous solid 

Upon incorporation of the additional velocity split, V s = V s + V? , described in Section [3.3.21 the 
resulting mass transport equation ([3]) for the solute species is 



dp s 
~d7 



7T 



div 



m s + 



P 



m 



p s div[v]. 



(15) 



In the hyperbolic limi t, where advection dominates, spatial oscillatio ns emerge in numerical solu- 
tions of this equation ([Brooks and Hug hel Il982l : iHughes et all Il987h . However, the form in which 
the equation is obtained is not in standard advection- diffusion form, and t herefore is not amenable 
to the application of standard stabilisation techniques ( Hughes et al.l . ll987l ). In part, this is because 
although the (near) incompressibility of the fluid phase is embedded in the balance of linear mo- 
mentum via the fluid stress, it has not yet been explicitly incorporated into the transport equations. 
This section proceeds to impose the fluid incompressibility condition and deduces implications for 
the solute mass transport equation, including a crucial simplification allowing for its straightforward 
numerical stabilisation. 

From Equation ([3]), the local form of the balance of mass for the fluid species (recalling that 
n f = 0) in the current configuration is 



dp* 



-div 



m 



/d 



iv v 



(16) 



In order to impose the incompressibility of the fluid, we first denote by Pq. . the initial value of 
the fluid reference concentration. Recall that the fluid concentration with respect to the reference 
configuration evolves in time; p$ = p (X,t). Therefore we can precisely, and non-trivially, define 

pIJ x ) 



pl(X,0) =:p f 0ini (X) 

= f & i (xo<p)J(X,t) 

p l {xOLp,t) 



jf*(X,t) 



-J(X,t) 



(17) 



pHx 



1 V t 



0(p ,t)j^lx,t). 



7 And not in terms of acceleration relative to fixed stars for e.g., as discussed in (|Truesdell and Noll 1 19651 . Page 
43). 

8 Where every observer has an implicit knowledge of the directionality of the field relative to a fixed frame, allowing it 
to transform objectively. Specifically, under a time-dependent rigid body motion imposed on the current configuration 
carrying x to x + = c(t) + Q(t)x, where c(t) £ R 3 and Q(t) £ SO(3), it is understood that the acceleration due to 
gravity in the transformed frame is g + — Q T g and is therefore frame-invariant. However, a + = c + 2Qv + Qx + Qa 
, and is therefore not frame-invariant. 
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In (|17p . J := det(.F) and J^ g := det(.F g ). The quantity p[ ni is denned by the right hand-sides 
of the first and second lines of (|17p . To follow the argument, consider, momentarily, a compressible 
fluid. If the current concentration, changes due to elastic deformation of the fluid and by 
transport, then pf ni as defined is not a physically-realized fluid concentration. It bears a purely 
mathematical relation to the current concentration, p f , since the latter quantity represents the 
effect of deformation of a tissue point as well as change in mass due to transport at that point. If 
the contribution due to mass change at a point is scaled out of p l the quotient is identical to the 
result of dividing p^. . by the deformation only. This is expressed in the relation between the right 
hand-sides of the second and third lines of (117p . The elastic component of fluid volume change 
m a pore is jf° := det(F c ), which appears in the third line of (|17p via the preceding arguments. 
Clearly then, for a fluid demonstrating near incompressibility intrinsically (i.e., the true density is 
nearly constant), we have J^ c ~ 1 as indicated. Equation (|17p therefore shows that for a nearly 
incompressible fluid occupying the pores of a tissue, if we further assume that the pore structure 
deforms as the solid collagenous skeleton, Pq(X,0) ps p f (x o ip,t). The fluid concentration as 
measured in the current configuration is approximately constant in space and time. This allows us 
to write, 



J, 



x o 



X 



0. 



(18) 



which is the hidden implication of our assumption of a homogeneous deformation, i.e., F is the 
deformation gradient of solid collagen and the pore spaces. This leads to = 00 We therefore 
proceed to treat our fluid mass transport at steady state. Rewriting the flux m f from Equation (fl~6j) 
as the product pV and using the result derived above, 







dpf 
at 

— div 



X 



pf v f 



(19) 



//div [v] 



Returning to (fl~5|) with this result, 



dff 
dt 




+ 7r s - div 



(20) 



Thus, using the incompressibility condition (|19p . we get the simplified form of the balance of 
mass for an arbitrary solute species, s, 



dff 

~dJ 



7r s - div 



rrr 



■ grad [p s ] p s mf • grad [p 



/l 



+ 



(21) 



Using the pushed- forward form of (|14p . this is now in standard advection-diffusion form, 



9 Which results in a very large pressure gradient driven flux due to incompressibility. 
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v v 

Advection term Additional, p s -dependent source term 



where D s is a positive semi-definite diffusivity, m/ j pf is the advective velocity, and 7r s is the 
volumetric source term. This form is well suited for sta bilisation schemes suc h as the streamline 
upwind Petrov-Galerkin (SUPG) method (see, for e.g., Irlughes et al. ( 19871 )). described briefly 
below, which limit spatial oscillations otherwise observed when the element Peclet number is large. 



3.3.5 Stabilisation of the simplified solute transport equation 

In weak form, the SUPG-stabilised method for Equation (|22p is 



h / d /° s f i 

w — h rn ■ grad 

n V dt 



+ 



L ( grad 


> 


■ D s grad 


P 












e=l 
n el 



mf 
t — j- ■ grad 
n c P 1 



'n, P i 



e=l 



W 



W 



dp s 



dt 



+ m,f ■ grad 



div 



D s grad 



w h ir s dQ+ w h h dT 



n el 



e=l 



+ T ~T ' grad 



n c P 



w 



tt s dn, 



dn 



(23) 



dn 



where quantities with the superscript h represent finite-dimensional approximations of infinite- 
dimensional field variables, r n is the Neumann boundary, and this equation introduces a numerical 
stabilisation p arameter r, which we have calculated from the L2 norms of element level matrices, 
as described in iTezduvar and Sathel hQ()A ). 



3.4 Nature of the sources 

There exists a large body of literature, ( Cowin and Hegedusl . 1976 ; Epstein and Maugin . 2000; 
Ambrosi and Mollical . 120021 ). that addresses growth in biological tissue mainly based upon a single 
species undergoing transport and production/annihilation. However, when chemistry is accounted 
for, it is apparent that growth depends on cascades of complex biochemical reactions involving sev- 
eral species, and additionally involves intimate coupling between mass t ransfer, biochemis t ry an d 
mechanics. An example of this chemo-mechanical coupling is described in |Proven Z anoetalJ &. 

The modelling approach followed in this work is to select appropriate functional forms of the 
source terms for collagen, IT, and the solutes, II s , that abstract the complexity of the biochemistry. 



15 




Figure 7: Engineered tendon constructs. See Calve et al.1 (|2004l ) for details 



In our earlier exposition ( Garikipati et al. . 20041 ) . we used simple first order chemical kinetics to 
define Il c . Other forms, which have been studied in th e literature, can be u sed: 

(i) Michaelis-Menten enzyme kinetics (see, for e.g., Sengers et al. ( 20041 )). which involves a two- 
step reaction with the collagen and solute production terms given by 



IT 



where p ce n is the concentration of fibroblasts, k. 



Pcell, 



it 



-IT 



(24) 



is the maximum value of the solute production 



reaction rate constant, and p s m is half the solute concentration corresp onding to k mKx . For detail s 



on the chemistry modelled by the Michaelis-Menten model, see, for e.g., Bromberg and Dill ( 20021 ). 

(ii) Strain energy- dependent sources that induce growth at a point when the energy density 
deviates from a reference val ue. An example of source term s of this form was originally proposed in 
the context of bone growth ( Harrigan and Hamilton . 19931 ). We are not aware of studies that have 
developed similar functional forms for soft tissue, and therefore have adapted this example from the 
bone growth literature, recognizing that this topic is in need of further study. Suitably weighted by 
a relative concentration ratio, and written for collagen, this source term has the form 



IT 



Po 

An, 



F ' 



(25) 



where ip-p is the mass-specific strain energy function, and ipp is a reference value of this strain energy 
density. Equation (|25p models collagen production when the strain energy density (weighted by a 
concentration ratio) at a point exceeds this reference value, and models annihilation otherwise. 



4 Numerical examples 

The theory presented in the preceding sections results in a system of non-linear, coupled partial 
differential eq uations. A fin it e element formulation em ploying a staggered scheme base d upon 
operator splits Armero ( 19991 ); Garikipati and Rao ( 200ll ) has been implemented in FEAP ( Taylor . 
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Figure 8: The finite element mesh used in the computations. 



19991 ) to solve the coupled problem. As an example, in the biphasic problem involving solid and fluid 



phases only, the basic solution scheme involves keeping the displacement field fixed while solving 
for the concentration fields using the mass transport equations. The resulting concentration fields 
are then fixed to solve the mechanics problem. This procedure is repeated until the resulting fields 
satisfy the differential equations within a specified numerical tolerance. 

The following examples aim to demonstrate the mathematical formulation and aspects of the 
coupled phenomena as the ti ssue grows. The m odel geometry, based on our engineered tendon 



constructs (see Figure [7J and I Calve et al.l (|2004h ). is a cylinder 12 mm in length and 1 mm in 



cross-sectional area. The corresponding finite element mesh using hexahedral elements, is shown in 
Figure El 

The following numerical examples involve solution of a common set of partial differential equa- 
tions. The constutive models, however, vary as we demonstrate the behaviour engendered by the 
many modelling assumptions discussed in the paper. The balance of linear momentum that we 
solve is ([6|) summed for 1 = c, f , with the constraint in ([7]) imposed. The absence of significant 
acceleration in the problems under consideration allows us to solve the balance of linear momentum 
quasi-statically. The fluid mass balance equation is solved in the current configuration, i.e. Q for 
i = f , but mass balance for the solid collagenous phase is solved in the reference configuration, i.e. 
dl]) for t = c. Mass balance for the solute is also solved in the current configuration, but using the 
stabilized scheme in weak form (|23p . The Backward Euler algorithm is used for all mass transport 
equations. The constitutive relation for the solid collagen follows ([8]). The constitutive relation for 
the fluid stress follows (El) with 



Kp { ) = U ( - 1 1 , (26) 




where K f is the fluid bulk modulus. The tissue is modelled as being fluid saturated in at t = 0, i.e. 
dSl) holds with pg sat = Po ini ■ However, the tissue is allowed to become unsaturated in for t > 
due to void formation. Then, the conditions set out in (llip apply. The chemical potential is then 
given by (|13p . The numerical examples that follow discuss further specialization of the constitutive 
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Parameter (Symbol) 


Value 


Units 


Chain density (N) 


7 x 1(P 




Temperature (9) 


310.6 


K 


Persistence length (^4) 


2.10 





Fully-stretched length (L) 


2.125 





Unit cell axes (a, b, c) 


1.95, 1.95, 2.43 





Bulk compressibility factors (7, (5) 


1000, 4.5 





Fluid bulk modulus (k/) 


1 


GPa 


Fluid mobility tensor (Dh = D f 5ij) 


1 x 10" 14 


s 


Fibroblast concentration (p ce ii) 


0.2 


kg.m -3 


Max. reaction rate (/c max = 5) 


5 


s" 1 


Max. solute concentration (p^J 


0.2 


kg.m -3 


Solute diffusivity (D s ) 


1 x 10~ 9 


m s 



Table 1: Material parameters used in the analysis. 



relations to other cases discussed in the preceding sections. The numerical values of parameters^! 
that have been used appear in Tabl e [p 



Non-linear projection methods (IS imp et all 11985) are used to treat the ne ar-incompressibility 



imposed by the fluid. Mixed methods, as described in lGarikipati and Raol (|200lh . are used for stress 
(and strain) gradient driven fluxes. 

The initial and boundary conditions have been chosen in order to model a few common me- 
chanical and chemical interventions on engineered tissue. However, we will not attempt detailed 
descriptions of experiments, choosing to focus instead on results that can be directly related to the 
models. A more detailed comparison with experiments is forthcoming in a separate communication. 



4.1 A multiphasic problem based on enzyme-kinetics 

This first example can be viewed as a model for localised, bolus delivery of regulatory chemicals 
to the tendon while accounting for mechanical (stress) effects. A single solute species 11 ! is consid- 
ered, denoted by s, and a uniform distribution of fibroblasts that are characterised by their cell 
concentration, p ce ii. Both, Fickean diffusion of the solute, and stress gradient driven fluid flow are 
incorporated in this illustration. We use Michaelis-Menten enzyme kinetics [Equation (|24|) ] to deter- 
mine the rates of solute consumption and collagen production as a function of solute concentration. 
This non-linear relationship for our choice of parameters is visualised in Figure EE Here, the fluid 
phase does not take part in reactions, and hence n f = 0. 

The tendon immersed in the bath is subjected to a constrictive radial load, such as would be 
imposed upon manipulating it with a set of tweezers, as depicted in Figure [TUl The maximum 
strain in the radial direction — experienced half-way through the height of the tendon — is 10%. The 
applied strain in the radial direction decreases linearly with distance from the central plane, and 
vanishes at the top and bottom surfaces of the tendon. 

I The mobility tensor reported in Table[T]is an order-of-magnitude estimate recalculated from lHan et all (j2000h to 
correspond to the mobility used in this paper. These authors reported a mean value of 0.927 x 10" 14 s, with a ranee of 

1.14 x 10~ 14 0.58 x 10 -14 s in terms of the mobility used here. Theirs is the mobility parallel to the fiber direction 

in Rabbit Achilles tendon. Our usage of it is as an isotropic mobility. Using anisotropic mobilities, or different values 
from the reported range changes the result quantitatively, but not qualitatively. 

II Here, we envision the solute to be a protein playing an essential role in growth by catalysing underlying biochemical 
reactions. An important example of this is a family of protein s, TGF/3, which is a multi-functional peptide that controls 
numerous functions of many cell types (|Alberts et all l2002h . 
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Figure 9: Variation of the collagen source term (kg.m 3 .s 1 ) with solute concentration (kg.m 3 ). 

N -M f 



The initial collagen concentration and the initial fluid concentration are both 500 kg.m -3 at 
every point in the tendon, and the fluid concentration in the bath is 500 kg.m" 3 . In addition, 
a solute-rich bulb of radius 0.15 mm is introduced with its centre on the axis of the tendon and 
situated 3 mm below the upper circular face of the tendon. The initial solute concentration is 
0.05 kg.m -3 at all other points in the tendon, and increases linearly with decreasing radius in this 
bulb to 1 kg.m -3 at its centre (see Figure [TTJ. The parameters used are listed in Table [H and are 
relevant to tendons. 

The aim of this example is to compare the influences upon solute transport from two mecha- 
nisms: Fluid stress gradient-driven transport, arising from the applied constrictive load, and solute 
concentration gradient-driven transport. These mechanisms have both been implicated in nutrient 
supply to cells in soft tissue. The results of this numerical example demonstrate that because the 
magnitude of the fluid mobility for stress gradient driven transport is orders of magnitude smaller 
than the diffusion coefficient for the solute through the fluid, there is relatively only a small stress 
gradient driven flux, and the transport of the solute is diffusion dominated. As a result, the solute 
diffuses locally, but displays no observable advection along the fluid. As the diffusion-driven solute 
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Figure 10: Constrictive load applied to tendon immersed in a bath. 



19 




Figure 11: The solute concentration (kg.m 3 ) initially. 




Figure 12: The collagen concentration (kg.m 3 ) at time t = 5 x 10 2 s. 

concentration in a region increases, the enzyme-kinetics model results in a small source term for 
collagen production, and we observe nominal growth. Figure [12] shows the collagen concentration 
at an early time, t = 5 x 10~ 2 s. 

This example incorporates all of the theory discussed in the paper. However, it is a valuable 
exercise in modelling to simplify the boundary value problem, and supress some of the coupled 
phenomena in order to gain a better understanding of some effcts. This is the approach followed in 
the next two numerical examples. The detailed transport and mechanics induced by the constrictive 
radial load are discussed first in Section 14.2. li 

4.2 Examples exploring the biphasic nature of porous soft tissue 

In these calculations, only two phases — fluid and collagen — are included for the mass transport and 
mechanics. The parameters used in the analysis are presented in Table [TJ 

4.2.1 The tendon under constriction 

In this example, the tendon immersed in a bath is subjected to the same constrictive radial load as in 
Section 14.11 Since that example demonstrated an insignificant amount of local collagen production 
over this time scale, we have simplified the problem by setting the source term II C = 0. The total 
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duration of the simulation is 10 s, and the radial strain is applied as a displacement boundary 
condition, increasing linearly from no strain initially to the maximum strain at time t = 1 s. Again, 
both the initial collagen concentration and the initial fluid concentration are 500 kg.m" 3 at every 
point in the tendon. This tendon is exposed to a bath where the fluid concentration is 500 kg.m -3 . 

While solving the balance of momentum for the biphasic problem of the solid collagen and a 
fluid phase, we currently treat the tissue as a single entity and employ a summation of Equation (|6j) 
over both species. Additionally, condition allows us to avoid constitutive prescription of the 
momentum transfer terms between solid collagen and fluid phases, q c and q f . This facilitates 
considerable simplification of the problem, but such a treatment requires additional assumptions on 
the detailed deformation of the constitutive phases of the tissue. An explicit assumption we have 
drawn on thus far is the equality of the deformation gradient of the solid collagen and pore spaces, 
allowing us to use the deformation gradient F, suitably decomposed to account for change in fluid 
concentration, to model the fluid stress. This assumption and its consequences have been discussed 
in Sections 12.11 12.21 12.2. 1} 13.21 13.2.11 and 13.3.41 Since the imposition of a common deformation 
gradient results in an upper bound for the effective stiffness of the tissue and magnitudes of the 
fluxes established, we refer to it as the upper bound model. This assumption plays a fundamental 
role in determining the fluid flux driven by the fluid stress gradient. 




Fluid flux fkQ/m"2/5) 



i 



-2.61 E-08 
-2.17E-08 
-1.74E-08 
-1.30E-08 
-8.70E-09 
-4.35E-09 
-5.59E-22 
4.35E-09 
8.70E-09 
1.30E-08 
1 .74E-08 
2.17E-08 
2.61 E-08 



Figure 13: Upper bound fluid flux (kg.m 2 .s 1 ) in the vertical direction at time t = 1 s. 

For this upper bound model, Figure [13] shows the fluid flux in the vertical direction at the final 
stage of the constriction phase of the simulation, i.e. at time t = 1 s. The flux values are positive 
above the central plane, forcing fluid upward, and negative below, forcing fluid fluid downward. 
This stress-gradient induced fluid flux results in a reference concentration distribution of the fluid 
that is higher near the top and bottom faces, as seen in Figure HH 

As a result, these regions would have seen a higher production of collagen, or preferential growth, 
in the presence of non-zero source terms. As discussed in Section [2. 1.1 1 the mass transport equations 
are solved in the current configuration, where physical boundary conditions can be set directly. 
The values reported in Figure Q3] are pulled back from the current configuration. The current 
concentrations do not change for this boundary value problem. 

Solving a problem of this nature in the reference configuration using p^ = const, as the boundary 
condition to represent immersion of the tendon in a fluid bath yields non-physical results, such as 
an unbounded flow. This occurs since the imposed strain gradient causes a stress gradient in the 
fluid that does not decay. The imposed boundary condition on Pq prevents a redistribution of 
concentration that would have provided an opposing, internal gradient of stress, which in turn 
would drive the flux to vanish. 
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Fluid Cone.. [kg/m A 3) 



4.89E+02 
4.80E+02 
4.82E+02 
4.84E+02 
4.86E+02 
4.88E+02 
4.90E+02 
4.92E+02 
4.94E+02 
4.96E+02 
4.98E+02 
5.00E+02 
5.00E+02 



Figure 14: Reference fluid concentration (kg.m 3 ) at time t = 1 s. 

The tendon is held fixed in the radial direction after the constriction phase. The applied stress 
sets up a pressure wave in the fluid travelling toward the top and bottom faces. As the fluid leaves 
these surfaces, we observe that the tendon relaxes. This is seen in Figure UH\ which plots the 
vertical displacement of the top face with time, showing a decrease in height of the tendon after the 
constriction phase. We keep the centre of the bottom face of the tendon fixed. 




Vertical displacement of trie top face with time ~ 



2 4 6 8 10 

Time (s) 



Figure 15: Relaxation of the top face of the tendon after the constriction phase. 

In order to define a range of the magnitude of fluid flux, we now introduce the lower bound 
model (on effective stiffness of the tissue and, consequently, the magnitude of the fluid flux). For this 
lower bound, we replace the earlier strain homogenisation requirement with a stress homogenisation 
requirement, viz. equating the hydrostatic stress of the solid phase and the fluid pressure in the 
current configuration: 



P 



3 L 



(27) 



where p i is the fluid pressure in the current configuration, tr[«] is the trace operator, and <x c = 
j^P c F c is the Cauchy stress of the solid. The Cauchy stress of an ideal fluid can be defined from 
its current pressure as a? — pi. Figure [16] reports the value of the vertical flux under the lower 
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bound modelling assumption, using boundary conditions identical to the previous calculation at 
time t = 1 s, the final stage of the constriction phase of the simulation. 




Time = 1.00E+00 s 



Figure 16: Lower bound fluid flux (kg.m 2 .s 1 ) in the vertical direction at time t = 1 s. 

The fluid flux values reported in Figures [TBI and [TBI (corresponding to the upper and lower bound 
modelling assumptions, respectively) are qualitatively similar, but differ by about three orders of 
magnitude. This wide range points to the importance of imposing the appropriate mechanical 
coupling model between interacting phases. Note, however, that we have computed bounds for the 
range of possible fluid flux values under the specified mechanical loading. Recall, furthermore, that 
the example in Section B~T1 used the upper bound model, and yet resulted in no discernible advective 
solute transport. This suggests strongly that, given the parameters in Tabled! convective transport 
of nutrients in tendons is dominated by diffusive transport. In future work, we will detail models 
that result in precise field values for the fluxes, which will replace the upper and lower bounds 
discussed here. 

This numerical example also points to the fact that a convenient measure of the strength of 
coupling between the mechanics and mass transport equations is the ratio of the variation in hy- 
drostatic stress of the fluid to that of the solid. In the lower bound case, where the fluid response 
is defined by Equation (f27|) . it is instructive to note that this ratio is unity. As a result, it is seen 
that the lower bound case exhibits significantly weaker coupling than the upper bound case. In 
the latter, variation in the common deformation gradient, 5F, causes instantaneous variation in 



5p i 0(k { 5F : F T ) and in i<5tr[<x c ] 0(k c SF : F T ), where k c is the bulk modulus of the solid. 



The strength of coupling between the equations plays a principal role in the rate of convergence 
of the solution, as observed in Table [21 where the residual norms of the equilibrium equation (and 
corresponding CPU times in seconds for an Intel(R)Xeon 3.4 GHz machine) are reported for the 
first 8 iterations of each of the two cases. Recall that the staggered scheme involves solution of the 
mechanics equation keeping the concentrations fixed, and the mass transport equation keeping the 
displacements fixed, in turn, until the solution converges. The table does not report the value of the 
residual norms arising from the solution of the mass transport equation for the fluid, which occurs 
after each reported solve of the of the mechanics equation. Although the initial mechanics residual 
norms in successive passes are decreasing linearly in both cases, the rapid decrease in this quantity 
in the weakly-coupled case ensures convergence in far fewer iterations than the strongly coupled 
case. Thus, the corresponding CPU times reported are also lower for the weakly coupled case. This 
is advantageous. In addition to being more physical, as argued at the beginning of Section 14.2.21 
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immediately below, the lower bound, weakly-coupled case makes it feasible to drive problems to 
longer, physiologically-relevant time-scales through the use of larger time steps. 

4.2.2 A swelling problem 



ic. (kq/m A 3) 

4.50E+02 
4.58E+02 
4.67E+02 
4.75E+02 
4.83E+02 
4.92E+02 
5.00E+02 
5.08E+02 
5.17E+02 
5.25E+02 
5.33E+02 
5.42E+02 
5.50E+02 

= 0.00E+00S 

Figure 17: The collagen initial concentration (kg.m -3 ). 



ic. (kq/m*3) 

5.85E+02 
5.89E+02 
5.92E+02 
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6.02E+02 
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6.09E+02 
6.12E+02 
6.15E+02 
6.19E+02 
6.22E+02 
6.26E+02 

1.80E+03S 

Figure 18: The collagen concentration (kg.m -3 ) after 1800 s. 

Motivated mainly by the recognition that the lower bound model for solid-fluid mechanical 
coupling ensures convergence to a self-consistent solution in just a few passes of the staggered 
solution scheme, we adopt this version of the coupling for our final problem. On this note we point 
out that solution of the individual balances of linear momentum equation for the solid collagenous 
and fluid phases with the momentum transfer terms [q c ,q^ in (0)] is a statement of momentum 
balance between them. There is reason to suppose, therefore, that equating the solid collagen and 
fluid stress, or some component of these tensors as done in the lower bound model, is a reasonable 
approximation to explicitly solving the balance of linear momentum for each phase, including the 
momentum transfers. In contrast, equating the deformation gradient of the solid collagen with 
deformation of the pore spaces subjects the fluid to a stress state also determined by this deformation 
gradient in the upper bound model. This approximation does not correspond to an underlying 
physical principle comparable to the satisfaction of individual balances of linear momentum for 
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Table 2: Mechanics equation residual norms and corresponding CPU times in seconds for the first 
8 passes of each of the two cases for a typical time increment, At = 0.1 s. 
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Figure 19: The stress (Pa) vs stretch curves before and after growth. 
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Figure 20: The volume of the tendon (m 3 ) evolving with time. Note the fluid transported-dominated 
regime until 25 s, followed by the longer reaction-dominated growth stage. 



solid collagen and fluid, with momentum transfers. It is therefore somewhat less motivated and 
more questionable. Clearly, a rigorous analysis or numerical comparisons of all three models: upper 
bound, lower bound and direct solution of individual solid-fluid momentum balances, must be carried 
out to conclusively demonstrate this. It is a possible topic for a future paper. 

In this example we will demonstrate the mechanical effects of growth due to collagen production. 
In the interest of focusing on this issue we assume that fibroblasts are available, and that the 
fluid phase bears the necessary nutrients for collagen production dissolved at a suitable, constant 
concentration. Collagen production is assumed to be governed by a first-order rate law. Newly- 
produced collagen has proteoglycan molecules bound to it, and they in turn bind water. We model 
this effect by associating a loss of nutrient-bearing free fluid with collagen production. A fluid sink 
II is introduced following first order kinetics, 



as m 



Carikipati et all (j2004h . Here k { is the reaction rate, taken to be 0.07 



(28) 

s _1 , and pQ. n . is the 

initial concentration of fluid. The collagen source is mathemaically equivalent to the fluid sink: 
II C = — II f . When p { Q > Po ini , collagen is produced. 

The boundary conditions in this example correspond to immersion of the tendon in a nutrient- 
rich bath. The initial collagen concentration is 500 kg.m -3 and the fluid concentration is 400 kg.m~ 3 
at every point in the tendon. When this tendon is exposed to a bath where the fluid concentration 
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is 410 kg.m -3 , i.e. p f (x, t) = 410 kg.m~ 3 Va; € dftf, nutrient-rich fluid is transported into the tissue, 
due to the pressure difference, induced by the concentration difference, between the fluid in the 
tendon and in the bath (fluid stress gradient-driven flux). Thereby, the nutrient concentration is 
elevated, leading to collagen production, fluid consumption and, eventually, growth due to additional 
collagen. 

Figure [T71 shows the initial collagen concentration in the tendon. After it has been immersed in 
the nutrient-rich bath for 1800 s, the tendon shows growth and the collagen concentration is higher 
as seen in Figure [18] On performing a uniaxial tension test on the tendon before and after growth, 
it is observed (Figure [T9|) that the grown tissue is stiffer and stronger due to its higher collagen 
concentration. Also note that there is a rapid, fluid transport-dominated swelling of the tendon 
between and 25 s following immersion in the fluid bath (Figure |20|) . This causes a small volume 
change of ~ 1.6%. In this transport-dominated regime the contribution to tendon growth from 
collagen production is small. However, the fluid-induced swelling saturates, and between 25 and 
1800 s the reaction producing collagen dominates the growth process, producing a further ~ 6.8% 
volume change. Noting that the range of collagen concentration in Figure [181 is 585—626 kg.m~ 3 , and 

that @ gives F &c = ( -£ 2 - ) 1, this portion of the volume change is quite clearly due to collagen 

production. The total volume change of 8.4% corresponds to changes in each linear dimension of 
the tendon by only ~ 2.7%, and is not discernible upon comparing Figures [T71 and [181 It is, however, 
manifest in Figure [20j 



5 Conclusion 

In this p aper, we have discussed a number of enhancements to our original growth formulation pre- 
sented in iGarikipati et al.l (j2004h . That formulation has served as a platform for posing a very wide 



range of questions on the biophysics of growth. Some issues, such as saturation, incompressibility of 
the fluid species and its influence upon the tissue response, and the roles of biochemical and strain 
energy-dependent source terms are specific to soft biological tissues. We note, however, that other 
issues are also applicable to a number of systems with a porous solid, transported fluid and reacting 
solutes. Included in these are issues of current versus reference configurations for mass transport, 
swelling, Fickean diffusion, fluid response in compression and tension, cavitation and the strength 
of solid-fluid coupling.. 

These issues have been resolved using arguments posed easily in the framework derived in 



Garikipati et al.l ([20041 ). The interactions engendered in the coupled reaction-transport-mechanics 
system are complex, as borne out by the numerical examples in Section SJ We are currently 
examining combinations of sources defined in Section 13. 4[ and aim to calibrate our choices from 
tendon growth experiments. The treatment of these issues has led to a formulation more suited to 
the biophysics of growing soft tissue, making progress toward our broader goal of applying it to the 
study of wound healing, pathological hypertrophy and atrophy, as well as a study of drug efficacy 
and interaction. 
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